Introduction to tmap

Introduction

This workshop will provide you with an introduction to making maps in R using the R package tmap. tmap is a relatively simple, but powerful, map making package. The code syntax will look somewhat familiar to ggplot2 users. Although ggplot2 can be used for map making, tmap is specifically designed for that purpose, and offers greater capabilities such as handling a wider variety of spatial data types and easy interactive mapping.

We will cover:

• Mapping polygons and points
• Mapping rasters
• Adding scale bars and a north arrow
• Adjusting legends
• Making a grid of small maps (facets)

Basic knowledge of R and spatial data is assumed for this workshop.

Resources:

tmap recently underwent a major upgrade with the release of version 4. There are a huge number of options for making and customizing maps, and this session will only provide a introduction to some of them. Here are some great resources where you can get more examples:

• The official tmap website, which includes many examples.
• Elegant and informative maps with tmap. An in-progress books written by the developers of tmap.
• The Geocomputation with R: 9 Making Maps in R book chapter. The entire books is available free online and is an excellent resource on map making and geospatial data in R in general.

Prerequisites

Apart from tmap, we will use the terra package for loading and manipulating spatial data. Both packages can be loaded using the following code

install.packages(c("terra", "tmap", "rnaturalearth"))

We can now load the packages:

library(terra)

terra 1.8.86

library(tmap)
library(rnaturalearth) #for adding country boundaries

Getting started with polygons

We are going to make a map using data for the Great Barrier Reef (GBR) Marine Park in Australia. First we will load some geospatial data we will need.

#load a geopackage with boundary of the Great Barrier Reef marine park
gbr_boundary <- vect("../data/gbr_MPA_boundary.gpkg")

#load a raster with sea surface temperature for the Great Barrier Reef
gbr_sst <- rast("../data/gbr_sst.tif")

Now we have some data loaded, lets jump straight into some map making. tmap is like the ggplot package: you have to tell it the layer you want to map, using tm_shape() and how you want to map that layer. Let’s try this with the data we have just loaded:

tm_shape(gbr_boundary) +
  tm_polygons()

We have a map! It shows the boundary of the Great Barrier Reef Marine Park. We used tm_polygons() because we want the polygon filled with a colour and a border drawn around it. You can also use tm_fill() or tm_borders() if you only want the polygon filled, or only the border. Here is a comparison:

map_fill <- tm_shape(gbr_boundary) +
  tm_fill() +
  tm_layout(panel.labels = "tm_fill()")

map_border <- tm_shape(gbr_boundary) +
  tm_borders() +
  tm_layout(panel.labels = "tm_borders()")

map_polygon <- tm_shape(gbr_boundary) +
  tm_polygons() +
  tm_layout(panel.labels = "tm_polygons()")

tmap_arrange(map_fill, map_border, map_polygon)

Interactive Map

Although we know our map is the GBR marine park boundary, we don’t have any other information on our map to get a sense of where it is. Sometimes we just want to see where our data are in the world! A quick and easy way to do that is using tmap’s interactive mode:

tmap_mode(mode = "view") #enter interactive mode

ℹ tmap mode set to "view".

  tm_shape(gbr_boundary) +
  tm_borders() 

tmap_mode("plot") #this takes us back into static map mode

ℹ tmap mode set to "plot".

I find this interactive mode very useful for exploring and zooming in on data, much like you might do in QGIS or ArcGIS software. You can use the interactive mode with any tmap map you make.

Adding raster data

Along with the boundary, we loaded some sea surface temperature data. We can just add this data onto the boundary we have already mapped:

tm_shape(gbr_boundary) +
  tm_borders() +
  tm_shape(gbr_sst) +
  tm_raster()

We mapped the sea surface temperature data using tm_raster() because it is raster data. But why can we only see a few parts of the GBR boundary now? This is because tmap stacks the layers in the order we list them. Since the raster was second, it was mapped on top of the boundary, obscuring most of it.

Let’s try again, changing only the order of the data:

tm_shape(gbr_sst) +
  tm_raster() +
tm_shape(gbr_boundary) +
  tm_borders() 

Better! We now see the boundary plotted over the temperature data.

Scales and legends

At the moment the temperature scale is in intervals of 2 degrees, starting at 18 and going up to 30. We might want to show it on a continuous scale and also change the colour palette to something more representative of low to high temperatures. We can control the way the values are displayed using the col.scale = argument to tm_raster().

Finding a suitable colour palette can be tricky, but I strongly recommend the cols4all package which has a great visual interface to help you choose a palette that not only looks good, but is colour-blind friendly. Use the command cols4all::c4a_gui() to access the interface. You will need to close it before you can continue running code.

tm_shape(gbr_sst) +
  tm_raster(col.scale = tm_scale_continuous(values = "ocean.thermal"),
            col.legend = tm_legend(title = "SST")) +
  tm_shape(gbr_boundary) +
  tm_borders()

I also changed the legend title using the col.legend = tm_legend(title = "SST"). The tm_legend(), like most tmap functions, has a lot of options (see ?tm_legend).

At the moment, our legend goes from high values at the bottom, to low values at the top. I generally prefer the opposite, and we can easily reverse the legend using reverse = TRUE:

tm_shape(gbr_sst) +
  tm_raster(col.scale = tm_scale_continuous(values = "ocean.thermal",
                                            ticks = seq(19, 30, by =2)),
            col.legend = tm_legend(title = "SST", 
                                   reverse = TRUE)) +
  tm_shape(gbr_boundary) +
  tm_borders()

When the legend is reversed, the marked values on the legend end up looking a bit odd because the lowest value marked on the scale is 22, even though temperatures go down to 19.4960003. To tidy the scale up, I manually set which values are shown on the legend using the ticks = argument to tm_scale_continuous()

Scale bar and compass

Our map looks pretty good now, but is missing a scale bar. You might also want to add a north arrow. There are lots of options for these (see ?tm_scale_bar and ?tm_compass), but we will just add the defaults, using just one extra line of code for each of these components:

tm_shape(gbr_sst) +
  tm_raster(col.scale = tm_scale_continuous(values = "ocean.thermal",
                                            ticks = seq(19, 30, by =2)),
            col.legend = tm_legend(title = "SST", 
                                   reverse = TRUE)) +
  tm_shape(gbr_boundary) +
  tm_borders() +
  tm_scalebar(position = c("bottom", "left")) +
  tm_compass(position = c("top", "right")) 

I have changed the default positions for the scale bar and compass to place them where there is space on the map.

Adding country boundaries

We often want to add country boundaries to a map. The rnaturalearth package provides various types of country boundaries at multiple resolutions. We will add polygons for the countries that fall within our map, using high resolution data from rnaturalearth, cropped to the area we are mapping:

countries_bordering <- rnaturalearth::ne_countries(scale = 'large', continent = "oceania", returnclass = "sv") |> #get data for Oceania as SpatVector (sv)
  crop(gbr_sst) #crop so that we only have countries within our mapping area

Now that we have the country polygons we need, we can add them to our map. We use the same process as before, using tm_shape() to specify the data, and we will use tm_polygons() to specify that we want polygons. The Natural Earth data has many columns (attributes) in it which contain information such as the population, G.D.P. and country names in many different languages. We will use the attribute “sovereignt” which is the sovereign state. There are only two states within the area we are mapping, so I will manually specify the colours to fill them, but you could also use a colour palette such as those found using cols4all::c4a_gui().

tm_shape(gbr_sst) +
  tm_raster(col.scale = tm_scale_continuous(values = "ocean.thermal",
                                            ticks = seq(19, 30, by =2)),
            col.legend = tm_legend(title = "SST", 
                                   reverse = TRUE)) +
  tm_shape(gbr_boundary) +
  tm_borders() +
  tm_shape(countries_bordering) +
  tm_polygons(fill = "sovereignt",
              fill.scale = tm_scale_categorical(values = c("tan", "forestgreen")),
              fill.legend = tm_legend(title = "")) +
  tm_scalebar(position = c("bottom", "left")) +
  tm_compass(position = c("top", "right")) 

Adding point data

We might also want to add some point data, such as survey locations, to our map. We will add some survey points from Crown-of-Thorns starfish surveys using the tm_dots() function to map the points, and remembering to use tm_shape() to tell tmap what data it is we want to use.

#Crown of Thorns survey data from https://data.aims.gov.au/data-download/ltmp-exports/nerp_ltmp_cots.csv

cots_data <- vect("../data/crown_of_thorns_survey.gpkg")

tm_shape(gbr_sst) +
  tm_raster(col.scale = tm_scale_continuous(values = "ocean.thermal",
                                            ticks = seq(19, 30, by =2)),
            col.legend = tm_legend(title = "SST", 
                                   reverse = TRUE)) +
  tm_shape(gbr_boundary) +
  tm_borders() +
  tm_shape(cots_data) + #this is the crown of thorns data
  tm_dots() +
  tm_shape(countries_bordering) +
  tm_polygons(fill = "sovereignt",
              fill.scale = tm_scale_categorical(values = c("tan", "forestgreen")),
              fill.legend = tm_legend(title = "")) +
  tm_scalebar(position = c("bottom", "left")) +
  tm_compass(position = c("top", "right")) 

We have mapped each survey point as a dot. Often, we will have information about how many species there are at a survey site. We can represent this using a colour scale, using the fill = argument to tm_dots(). Another option is to size the dots according to the number using size =.

tm_shape(gbr_sst) +
  tm_raster(col.scale = tm_scale_continuous(values = "ocean.thermal",
                                            ticks = seq(19, 30, by =2)),
            col.legend = tm_legend(title = "SST", 
                                   reverse = TRUE)) +
  tm_shape(gbr_boundary) +
  tm_borders() +
  tm_shape(cots_data) + #this is the crown of thorns data
  tm_dots(fill = "MEAN_COTS",
          fill.legend = tm_legend(title = "Mean No. COTS")) +
  tm_shape(countries_bordering) +
  tm_polygons(fill = "sovereignt",
              fill.scale = tm_scale_categorical(values = c("tan", "forestgreen")),
              fill.legend = tm_legend(title = "")) +
  tm_scalebar(position = c("bottom", "left")) +
  tm_compass(position = c("top", "right")) 

Facets: Many maps

Often we want to show a grid of maps that all look basically the same, but with one variable different in each one. For example, we might want to show the temperature for several months or years so we can visually compare the maps.

For raster data, we often store data for the same area in a multi-layer raster. As an example, we will load a multi-layer raster with monthly sea surface temperature for 12 months.

monthly_sst <- rast("../data/gbr_monthly_sst.tif")

monthly_sst

class       : SpatRaster 
size        : 301, 261, 12  (nrow, ncol, nlyr)
resolution  : 0.05, 0.05  (x, y)
extent      : 142, 155.05, -25, -9.95  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : gbr_monthly_sst.tif 
names       :  2022_06,  2022_07,  2022_08,  2022_09,  2022_10,  2022_11, ... 
min values  : 19.49600, 15.93677, 17.18129, 20.42967, 22.46194, 23.90833, ... 
max values  : 29.81233, 29.54355, 29.25194, 29.71567, 30.24129, 30.73400, ... 

We can see there are 12 layers (nlyr) and the names are the year and month of the data. What happens if we take the code we used above for mapping and just replace the single raster layer with our multi-layer raster?

tm_shape(monthly_sst) +
  tm_raster(col.scale = tm_scale_continuous(values = "ocean.thermal"),
            col.legend = tm_legend(title = "SST", 
                                   reverse = TRUE)) +
  tm_shape(gbr_boundary) +
  tm_borders() +
  tm_shape(countries_bordering) +
  tm_polygons(fill = "sovereignt",
              fill.scale = tm_scale_categorical(values = c("tan", "forestgreen")),
              fill.legend = tm_legend(title = ""))

[plot mode] fit legend/component: Some legend items or map compoments do not
fit well, and are therefore rescaled.
ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

Impressive! tmap automatically makes a grid of small maps (facets). I removed the scale bar and north arrow otherwise the maps would be a bit crowded, but you can leave them in if you like. I also removed the parts of the code that specified the limits and ticks since these will be different now we have more data.

Having a scale for each map is not ideal because would like to have all the maps on the same scale so we can compare the ocean temperatures across months. We can make it so that we only have one scale using col.free = FALSE in tm_raster().

tm_shape(monthly_sst) +
  tm_raster(col.scale = tm_scale_continuous(values = "ocean.thermal"),
            col.legend = tm_legend(title = "SST", 
                                   reverse = TRUE), 
            col.free = FALSE) +
  tm_shape(gbr_boundary) +
  tm_borders() +
  tm_shape(countries_bordering) +
  tm_polygons(fill = "sovereignt",
              fill.scale = tm_scale_categorical(values = c("tan", "forestgreen")),
              fill.legend = tm_legend(title = ""))

We have a single scale for all maps now and we can see the ocean gets warmer in the summer months (November to April - remember this is the Southern hemisphere!), and cooler in the winter. The numbers on our scale have gone a bit weird again though. To fix that, we can set the ticks = argument manually, but we will need to know the minimum and maximum temperature values.

#minimum temperature
min_temp <- global(monthly_sst, min, na.rm = TRUE) |> min()

#maximum temperature
max_temp <- global(monthly_sst, max, na.rm = TRUE) |> max()

#create sequence of numbers for the scale
scale_numbering <- seq(round(min_temp), round(max_temp), by = 2)

We can also create nicer month and year labels for the maps.

month_labels <- c(month.abb[6:12], month.abb[1:5])

month_year_labels <- paste(month_labels, c(rep(2022, 7), rep(2023, 5)), sep = " ")

We can add these labels using panel.labels = argument in the tm_layout() function:

tm_shape(monthly_sst) +
  tm_raster(col.scale = tm_scale_continuous(values = "ocean.thermal",
                                            ticks = scale_numbering),
            col.legend = tm_legend(title = "SST", 
                                   reverse = TRUE), 
            col.free = FALSE) +
  tm_shape(gbr_boundary) +
  tm_borders() +
  tm_shape(countries_bordering) +
  tm_polygons(fill = "sovereignt",
              fill.scale = tm_scale_categorical(values = c("tan", "forestgreen")),
              fill.legend = tm_legend(title = "")) +
  tm_layout(panel.labels = month_year_labels)